Assignment 4

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.2     ✓ purrr   0.3.4
✓ tibble  3.0.4     ✓ dplyr   1.0.2
✓ tidyr   1.1.2     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.0
── Conflicts ───────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(lubridate)

Attaching package: ‘lubridate’

The following objects are masked from ‘package:base’:

    date, intersect, setdiff, union

library(modelr)
library(broom)

Attaching package: ‘broom’

The following object is masked from ‘package:modelr’:

    bootstrap

Modeller

Leser inn data

arblos <- read_csv("data/al9914m.csv")

── Column specification ───────────────────────────────────────────────────
cols(
  knr = col_character(),
  knavn = col_character(),
  aar = col_double(),
  mnd = col_double(),
  al_Menn = col_double(),
  al_Kvinner = col_double(),
  alp_Menn = col_double(),
  alp_Kvinner = col_double(),
  alp_15_74 = col_double(),
  alp_15_29 = col_double(),
  alp_30_74 = col_double()
)

bef <- read_csv("data/bef9914MK.csv")

── Column specification ───────────────────────────────────────────────────
cols(
  knr = col_character(),
  knavn = col_character(),
  aar = col_double(),
  bef_K_0_14 = col_double(),
  bef_K_15_29 = col_double(),
  bef_K_30_74 = col_double(),
  bef_K_75_105 = col_double(),
  bef_M_0_14 = col_double(),
  bef_M_15_29 = col_double(),
  bef_M_30_74 = col_double(),
  bef_M_75_105 = col_double(),
  bef_MK_0_14 = col_double(),
  bef_MK_15_29 = col_double(),
  bef_MK_30_74 = col_double(),
  bef_MK_75_105 = col_double()
)

Modeller med data fra bef (befolkning)

dim(arblos)

[1] 77330    11

dim(bef)

[1] 6688   15

names(bef)

 [1] "knr"           "knavn"         "aar"           "bef_K_0_14"   
 [5] "bef_K_15_29"   "bef_K_30_74"   "bef_K_75_105"  "bef_M_0_14"   
 [9] "bef_M_15_29"   "bef_M_30_74"   "bef_M_75_105"  "bef_MK_0_14"  
[13] "bef_MK_15_29"  "bef_MK_30_74"  "bef_MK_75_105"

names(arblos)

 [1] "knr"         "knavn"       "aar"         "mnd"         "al_Menn"    
 [6] "al_Kvinner"  "alp_Menn"    "alp_Kvinner" "alp_15_74"   "alp_15_29"  
[11] "alp_30_74"  

Arbeidsledighetsprosenten blir beregnet som: arbl% = antall arb. ledige/arbeidsstyrken. Arbeidsstyrken er her dem man anser egnet for arbeid dvs. uføretrygdete etc. er trukket ut. Arbeidsstyren i en aldesrkategori er derfor langt mindre enn befolkningen i den tilsvarende alderskategorien.

Dessverre kjenne vi ikke arbeidsstyrken, men vi kan beregne den vha.: arbeidsstyrken = antall arb. ledige/arb. ledighetsprosent.

arblos <- arblos %>% 
    mutate(
        wf_K = (al_Kvinner/alp_Kvinner)*100,
        wf_M = (al_Menn/alp_Menn)*100,
        wf_KM = wf_K + wf_M
    )

arblos

Arbeidsstyrken

Når befolkningen øker vil også arbeidsstyrken øke. Det er derfor mer naturlig å se på arbeidsstyrken relativt til delen av befolkningen som er i yrkesaktiv alder (15-74 år her).

names(bef)

 [1] "knr"           "knavn"         "aar"           "bef_K_0_14"   
 [5] "bef_K_15_29"   "bef_K_30_74"   "bef_K_75_105"  "bef_M_0_14"   
 [9] "bef_M_15_29"   "bef_M_30_74"   "bef_M_75_105"  "bef_MK_0_14"  
[13] "bef_MK_15_29"  "bef_MK_30_74"  "bef_MK_75_105"

Vi skal starte med å lage et datasett med arbeidsstyrken (wf) for hele landet samlet, men fordelt på de tre kategorien kvinner, menn og kvinner + menn.

Bruk data for januar hvert år til å beregne wf på landsbasis

# årlige data landet samlet
wf <- arblos %>% 
    filter(mnd == 1) %>% 
    group_by(aar) %>% 
    summarise(
        wf_K = sum(wf_K, na.rm = TRUE),
        wf_M = sum(wf_M, na.rm = TRUE),
        wf_KM = wf_K + wf_M
    )

`summarise()` ungrouping output (override with `.groups` argument)

dim(wf)

[1] 16  4

wf

dim(wf)

[1] 16  4

names(wf)

[1] "aar"   "wf_K"  "wf_M"  "wf_KM"

Summer de ulike årskategoriene for de to kjønnene og menn+kvinner for å finne total befolkning de ulike årene. Bruk mutate til å lage de nye variablene.

 bef <- bef %>% 
  mutate(
    bef_K = bef_K_0_14 + bef_K_15_29 + bef_K_30_74 + bef_K_75_105,
    bef_M = bef_M_0_14 + bef_M_15_29 + bef_M_30_74 + bef_M_75_105,
    bef_KM = bef_K + bef_M)

dim(bef)

[1] 6688   18

names(bef)

 [1] "knr"           "knavn"         "aar"           "bef_K_0_14"   
 [5] "bef_K_15_29"   "bef_K_30_74"   "bef_K_75_105"  "bef_M_0_14"   
 [9] "bef_M_15_29"   "bef_M_30_74"   "bef_M_75_105"  "bef_MK_0_14"  
[13] "bef_MK_15_29"  "bef_MK_30_74"  "bef_MK_75_105" "bef_K"        
[17] "bef_M"         "bef_KM"       

Legg befolkningsdata varaiablene bef_K, bef_M, bef_KM til wf. Husk at de må aggregeres for hele landet (group_by() og så summarise() før de «joines»). Bruk tilslutt mutate() make_date() for å lage en ny variabel år som er en date, dvs aar + month=1L + day=1L.

wf_hjelp <- bef %>%
  group_by(aar) %>%
  summarise(bef_K = sum(bef_K),
            bef_M = sum(bef_M),
            bef_KM = sum(bef_KM))

`summarise()` ungrouping output (override with `.groups` argument)

wf <- left_join(wf, wf_hjelp)

Joining, by = "aar"

wf <- wf %>%
  select(aar, bef_K, bef_M, bef_KM, wf_K, wf_M, wf_KM) %>%
  mutate(år = make_date(aar, month = 1L, day = 1L))

rm(wf_hjelp)

Da skal du få

wf

Vi vil nå se på arbeidsstyrke relativt til befolkning på landsbasis, dvs. wf_K/bef_K etc. Plot dataen vha. geom_line() for de tre kategoriene. Bruk år som x-variabel.

wf %>% 
    ggplot(wf, mapping = aes(x = år)) +
  geom_line(mapping = aes(y = wf_K/bef_K), color = "red") +
  geom_line(mapping = aes(y = wf_M/bef_M, color = "blue"), color = "darkblue") +
  geom_line(mapping = aes(y = wf_KM/bef_KM), color = "black") +
  ggtitle("Arbeidsstyrke relativt til befolkning. Menn blå, kvinner rød og totalt svart.") +
  theme(legend.position = "none")

Hvordan kan pukkelen rett før 2010 forklares?

Det kan skyldes at verdensøkonomien gikk bra, og at det i 2008/2009 var en finanskrise i verden som påvirket økonomien og næringene.

names(arblos)

 [1] "knr"         "knavn"       "aar"         "mnd"         "al_Menn"    
 [6] "al_Kvinner"  "alp_Menn"    "alp_Kvinner" "alp_15_74"   "alp_15_29"  
[11] "alp_30_74"   "wf_K"        "wf_M"        "wf_KM"      

Vi vil nå generere dat for arbeidsstyrken på fylkesbasis. Husk at de to første sifferene i knr angir fylket en kommune ligger i. Bruk dataene i arblos til å finne arbeidsstyrken på fylkesbasis (wf_f). Start med å bruke mutate() til å lage en ny variabel fylke. Grupper så og finn wf_K_f, wf_M_f og wf_KM_f vha. summarise(). Lag til slutt en ny variabel år som en date. Velg til slutt de relevante variabelen vha. select().

# årlige data per fylke
wf_f <- arblos %>% 
    mutate(
        fylke = substr(knr, start = 1, stop = 2)
    ) %>% 
    group_by(aar, mnd, fylke) %>% 
    summarise(
        wf_K_f = sum(wf_K, na.rm = TRUE),
        wf_M_f = sum(wf_M, na.rm = TRUE),
        wf_KM_f = wf_K_f + wf_M_f
    ) %>% 
    mutate(år = make_date(aar, mnd, day=1L)) %>% 
    select(aar, mnd, år, fylke, wf_K_f, wf_M_f, wf_KM_f)

`summarise()` regrouping output by 'aar', 'mnd' (override with `.groups` argument)

Vi skal nå ha en tibble wf_f som ser slik ut:

print(wf_f, n = 5)

Lag også en ny tibble bef_f fra bef som inneholder befolkningen i hvert fylke hvert år.

bef_f <- bef %>% 
      mutate(
        fylke = substr(knr, start = 1, stop = 2)
    ) %>% 
    group_by(aar, fylke) %>% 
    summarise(
        bef_K_f = sum(bef_K, na.rm = TRUE),
        bef_M_f = sum(bef_M, na.rm = TRUE),
        bef_KM_f = bef_K_f + bef_M_f
    ) %>% 
    select(aar, fylke, bef_K_f, bef_M_f, bef_KM_f)

`summarise()` regrouping output by 'aar' (override with `.groups` argument)

bef_f

Lag til slutt tibble-en wf_f_bef som innholde arbeidsstyrke (wf-f) og befolkning (bef-f) på fylkesnivå for hvert år.

wf_f_bef <- wf_f %>% 
  right_join(bef_f) %>%
  arrange(fylke, aar, mnd)

Joining, by = c("aar", "fylke")

Da skal wf_f_bef se slik ut:

print(arrange(wf_f_bef, fylke, aar, mnd))

Plot nå arbeidsstyrke relativt til befolkning vha. geom_line()

# kvinner og menn
wf_f_bef %>%
  ggplot() +
  geom_line(wf_f_bef, mapping = aes(x = år,
                                    y = wf_KM_f/bef_KM_f,
                                    color = fylke)) +
  theme(legend.position = "bottom")

Lage regioner

Alle fylkene blir litt rotete så vi definerer istedet seks regioner vha. case_when() og lager en ny tibble wf_r fra wf_f_bef.

# region
wf_r <- wf_f_bef %>% 
    mutate(
         dato = ymd(paste(aar, mnd, "01", sep = "-")),
        region = case_when(
            as.numeric(fylke) == 3 ~ "Oslo",
            as.numeric(fylke) %in% c(1:2, 4:8) ~ "Østlandet",
            as.numeric(fylke) %in% c(9, 10) ~ "Sørlandet",
            as.numeric(fylke) %in% c(11, 12, 14, 15) ~ "Vestlandet",
            as.numeric(fylke) %in% c(16, 17) ~ "Trøndelag",
            as.numeric(fylke) %in% c(18, 19, 20) ~ "Nord-Norge"
        )
    ) %>% 
    group_by(dato, region) %>% 
    summarise(
        wf_K_r = sum(wf_K_f, na.rm = TRUE),
        wf_M_r = sum(wf_M_f, na.rm = TRUE),
        wf_KM_r = wf_K_r + wf_M_r,
        bef_K_r = sum(bef_K_f, na.rm = TRUE),
        bef_M_r = sum(bef_M_f, na.rm = TRUE),
        bef_KM_r = bef_K_r + bef_M_r
    ) %>% 
    select(dato, region, wf_K_r, wf_M_r, wf_KM_r, bef_K_r, bef_M_r, bef_KM_r)

`summarise()` regrouping output by 'dato' (override with `.groups` argument)

Plot nå for regionene wf_KM_r/bef_KM_r, både vha. geom_line() og geom_smooth(). La farge vise regionene. Sett denne i ggplot() slik at det gjelder for både geom_line() og geom_smooth(). Sett i tillegg alpha = 0.5 for geom_line() og se = FALSE for geom_smooth(). Legg til theme(legend.position = “bottom”) til slutt for å få legend under plottet.

# kvinner og menn
wf_r %>% 
  ggplot(wf_r, mapping = aes(x = dato, 
                             y = wf_KM_r/bef_KM_r,
                             color = region)) +
  geom_line(alpha = 0.5) +
  geom_smooth(se = FALSE) +
  theme(legend.position = "bottom") +
  ggtitle("Arbeidsstyrkens andel av befolkningen for regionene.")

Lag tilsvarende plot for kvinner.

# kvinner
wf_r %>% 
  ggplot(wf_r, mapping = aes(x = dato, 
                             y = wf_K_r/bef_K_r,
                             color = region)) +
  geom_line(alpha = 0.5) +
  geom_smooth(se = FALSE) +
  theme(legend.position = "bottom") +
  ggtitle("Arbeidsstyrkens andel av befolkningen for kvinner per region.")

Lag tilsvarende plot for menn

# menn
wf_r %>% 
  ggplot(wf_r, mapping = aes(x = dato, 
                             y = wf_M_r/bef_M_r,
                             color = region)) +
  geom_line(alpha = 0.5) +
  geom_smooth(se = FALSE) +
  theme(legend.position = "bottom") +
  ggtitle("Arbeidsstyrkens andel av befolkningen for menn per region.")

Forlar kort den generelle utviklingen i arbeidsstyrken. Hva er det som «redder oss»?

For kvinnene kan en se at det er en jevnt økende trend at flere kvinner fra arbeidsstyrken jobber. Det kan skyldes flere grunner, men en av dem kan være kampen for likestilling og at kvinnene skal aktivt delta i jobbmarkedene og ikke slik de gjorde før. En kan se at de har en liten utflating av kurven i 2007/2008 og den faller litt etter 2011. Den første utflatingen kan skyldes finanskrisen og det som førte med der, mens i 2011 kan det skyldes andre ting. En kan se at hos menn synker kurven i alle regionene, i noen regioner mer enn andre. Litt fargesvak, men tror det er Oslo som redder oss. Eventuelt Vestlandet.

nest() arblos

Da skal vi jobbe direkte med arbeidsløshet og lage lineære modeller for hver av de 418 kommunene. Modellen vi skal lage er på ingen måte perfekt. Vi er interessert i selve teknikken med å organisere dataene og kjøre modeller på mange subsett av dataene (her for hver kommune).

Vi vil se på en modell der vi forklarer arbeidsstyrken i en kommune vha. ungdomsledighet og ledighet blant litt eldre arbeidstakere (30-74 år). En hypotese er at vi vil se en negativ sammenheng mellom ungdomsledighet og arbeidsstyrken. De unge flytter hvis det ikke er jobb.

Vi starter med å gruppere på kommune og nest-e dataene.

arblos_by_knr <- arblos %>% 
    group_by(knr, knavn) %>% 
    nest()

print(arblos_by_knr, n = 4)

arblos_by_knr$data[[1]]

Vi har nå en tibble med data for hver kommune inne i tibble-en arblos_by_knr.

Skriv en funksjon som kjører den lineære modellen wf_KM ~ alp_15_29 + alp_30_74 på en input dataframe a_df. Kall funksjonen mod1 (i magel på noe bedre navn).

mod1 <- function(a_df){
  lm(wf_KM ~ alp_15_29 + alp_30_74, data = a_df)}

KJøre så modellen vha. funksjonen mod1 på data i arblos_by_knr, og lag en list-column i arblos_by_knr som inneholder modellen. Kjør også tidy og glance fra broom på modellene for å få hhv. koeffisienter og ulike summary av modellen. Lagre resultatene i hhv. mod1_arblos, mod1_arblos_coef og mod1_arblos_sum (i siste er sum forkortelse for summary)

arblos_by_knr <- arblos_by_knr %>% 
    mutate(mod1_arblos = map(data, .f = mod1)) %>%
    mutate(mod1_arblos_coef = map(mod1_arblos, .f = tidy)) %>% 
    mutate(mod1_arblos_sum = map(mod1_arblos, .f = glance))
arblos_by_knr[arblos_by_knr$knr == "0441",][[3]]

[[1]]

NA

arblos_by_knr skal nå se ut slik:

arblos_by_knr

Kommunestørrelse

Vi lager oss så en ny kategori variabel for kommunestørrelse.

kom_str <- bef %>% 
    filter(aar == 2014) %>% 
    mutate(
        k_str = case_when(
            bef_KM <= 2500 ~ "Svært liten",
            bef_KM > 2500 & bef_KM <= 6000 ~ "Liten",
            bef_KM > 6000 & bef_KM <= 20000 ~ "Middels",
            bef_KM > 20000 & bef_KM <= 60000 ~ "Stor",
            bef_KM > 60000 & bef_KM <= 1000000 ~ "Svært stor"
        )
    ) %>% 
    select(knr, knavn, bef_KM, k_str)
head(kom_str)

Pakker ut og henter model karakteristika

Bruk så unnest() til å pakke ut mod1_arblos_coef og mod1_arblos_coef (husk at de to må stå i en c() og ha anførselstegn). Plukk ut variablene knr, knavn, term, estimate, std.error, p.value…9, adj.r.squared og legg resultatet i mod_arbl_re.

#modell arbeidsløse resultat
mod_arbl_re <- arblos_by_knr %>% 
  unnest(c("mod1_arblos_coef", "mod1_arblos_sum"), names_repair = "universal") %>% 
  select(knr, knavn, term, estimate, std.error, p.value...9, adj.r.squared)

New names:
* statistic -> statistic...8
* p.value -> p.value...9
* statistic -> statistic...13
* p.value -> p.value...14

Du skal da ha noe som ser slik ut:

print(mod_arbl_re, n = 10)

Slå sammen kom_str og mod_arbl_re vha. left_join(). Kall resultatet for kom_str_mod.

kom_str_mod <- mod_arbl_re %>% 
    left_join(kom_str, by = c("knr", "knavn")) 

Da skal du ha noe som ser slik ut:

print(kom_str_mod, n = 5)

Vi plotter koeffisientene som er signifikante og lar farge vise kommune størrelse. Tar vekk ekstreme estimat.

kom_str_mod %>% 
    filter(term == "alp_15_29") %>% 
    filter(estimate > -500 & estimate < 500 ) %>%
    filter(p.value...9 < 0.05) %>% 
    ggplot(mapping = aes(
        x = estimate,
        y = p.value...9,
        colour = k_str
        )
        ) +
    geom_point()

Lag tilsvarende plot for alp_30_74.

kom_str_mod %>% 
    filter(term == "alp_30_74") %>% 
    filter(estimate > -500 & estimate < 500 ) %>%
    filter(p.value...9 < 0.05) %>% 
    ggplot(mapping = aes(
        x = estimate,
        y = p.value...9,
        colour = k_str
        )
        ) +
    geom_point()

Hvor mange har vi?

kom_str_mod %>% 
    filter(term == "alp_15_29") %>% 
    group_by(k_str) %>% 
    summarise(n = n())

`summarise()` ungrouping output (override with `.groups` argument)

kom_str_mod %>% 
    filter(term == "alp_15_29") %>% 
    filter(estimate > - 500 & estimate < 500 ) %>%
    filter(p.value...9 < 0.05) %>% 
    group_by(k_str) %>% 
    summarise(n = n())

`summarise()` ungrouping output (override with `.groups` argument)

kom_str_mod %>% 
    filter(term == "alp_30_74") %>% 
    filter(estimate > - 500 & estimate < 500 ) %>%
    filter(p.value...9 < 0.05) %>% 
    group_by(k_str) %>% 
    summarise(n = n())

`summarise()` ungrouping output (override with `.groups` argument)

kom_str_mod %>% 
    filter(term == "alp_15_29") %>% 
    filter(p.value...9 < 0.05) %>% 
    group_by(k_str) %>% 
    summarise(mean_15_29 = mean(estimate))

`summarise()` ungrouping output (override with `.groups` argument)

kom_str_mod %>% 
    filter(term == "alp_30_74") %>% 
    filter(p.value...9 < 0.05) %>% 
    group_by(k_str) %>% 
    summarise(mean_30_74 = mean(estimate))

`summarise()` ungrouping output (override with `.groups` argument)

kom_str_mod %>% 
    filter(term == "alp_15_29") %>% 
    filter(p.value...9 < 0.05) %>% 
    group_by(k_str) %>% 
    summarise(n = n())

`summarise()` ungrouping output (override with `.groups` argument)

kom_str_mod %>% 
    filter(term == "alp_30_74") %>% 
    filter(p.value...9 < 0.05) %>% 
    group_by(k_str) %>% 
    summarise(n = n())

`summarise()` ungrouping output (override with `.groups` argument)

I litt over 50% av kommunene, hovedsaklig de små, ser modellen ut til virke. Kanskje noe å bygge videre på.

#siste